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In this paper we show how the transposition, the basic operation of the permutation group, can 
be taken into account in a diffusion process of identical particles. Whereas in an earlier approach 
fH ^ the method was applied to systems in which the potential is invariant under interchanging the 

O ' Cartesian components of the particle coordinates, this condition on the potential is avoided here. 

O ' In general, the potential introduces a switching of the boundary conditions of the walkers. [These 

transitions modelled by a continuous-time Markov chain generate sample paths for the propagator 
as a Feynman-Kac functional . A few examples , including harmonic fermions with an anharmonic 
I interaction, and the ground-state energy of ortho-helium are studied to elucidate the theoretical 

discussion and to illustrate the feasibility of a sign-problem-free implementation scheme for the 
recently developed many-body diffusion approach.] 
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I. INTRODUCTION 



In a previous paper Q, the present authors have shown that the many-body diffusion algorithm (MB DA) allows 
for a sign-problem-free simulation of excited antisymmetric states of interacting harmonic oscillators. The basic 
idea behind this approach has been to split the Euclidean-time propagator into a sum of independent propagators 
which remain positive on the appropriate state space |^-^ . Each of the elementary propagators could be individually 
sampled as a diffusion process for distinguishable particles on a state space with absorbing or reflecting boundaries. 
In the present paper, the many-body diffusion formalism and its implementation are generalized to allow for the 
' construction of interdependent diffusion processes. In this approach the propagator that governs the Euclidean time 
T-H , evolution is split into two parts: the kinetic part is used to describe the evolution of the system with well-defined 
■"si" ■ boundary conditions, the potential is used for branching and killing and rules transitions from one set of boundary 
^5 ' conditions to another set. The evolution is given by a Brownian motion or a random walk, but the boundary conditions 
can change during the simulation according to a continuous-time Markov process. The transition rates of this process 
partly or completely derive from the potential which determines also the killing or branching process. 

This principle is illustrated with some test models. A limited number of degrees of freedom is considered to avoid 
large computation times and to allow for a visual control on the probability densities that are generated. The presented 
algorithm does not integrate smoothly with the previous one for harmonic fermions, in the sense that for this algorithm 
'"^ • the interchange of two identical particles is fully under control but in combination with the even permutations it is 
C5 I far from trivial if more than two particles of parallel spin are present. 

O ■ In the next section, the theoretical background of the method is discussed. In Sec. Ill, details on the diffusion 
in a reduced state space are provided. In Sec. IV, a sign-problem-free implementation is developed and applied to 
two toy models to illustrate the concepts. Subsequently, anisotropic harmonic fermions are studied avoiding normal 
coordinates in order to elucidate the role of the potential for the transitions between different sets of boundary 
conditions. Finally, the MBDA is employed to calculate the ground-state energy of ortho-helium. Some concluding 



X 

■ remarks are presented in Sec. V. 
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II. A SIGN-PROBLEM-FREE APPROACH 



In this section, we describe the construction of a Feynman-Kac functional equipped with a process to generate the 
sample paths over which the exponential containing the interaction has to be averaged. This construction is based 
on the many body diffusion approach introduced in |^] and illustrated in 0|. The procedure is sign-problem- free in 
the strict sense, i.e. expectation values in the Feynman-Kac functional are obtained using transition probabilities in 
a Wiener-Poisson space. A mathematical account on the foundations of these compound processes can be found in 
. The Wiener space is the support for the diffusion whereas a multivariate Poisson process counts the switches of 
the boundary condition as will be explained below. 



A. Diffusion on a Domain 



Let r be a 3A^-dimensional point in the configuration space IR. of N distinguishable particles and assume that the 
propagator for this system is given by: 



KD[rf,T;ri,0] = ^ (V'fc(r,),exp(-ri?) '(/'fe(r/)) , 



(1) 



where H describes the evolution in Euclidean time of the N particles. The subscript D denotes that at this stage 
the particles are still considered to be distinguishable. This propagator can be easily written as a Feynman-Kac 
functional over all paths starting in Ti and ending in r^^ a Euclidean time lapse r later. The paths are constructed 
by a 3iV-dimensional Brownian motion {i? (r) ; r > 0} with variance ~ hr/rn 



(2) 



All integration paths in (||) start in as indicated by the averaging index £V. and satisfy the condition i? (r) — rf 
as denoted by the indicator I ^-^(^^-^^-^ Writing the propagator as an average over a Brownian motion is equivalent 
to giving the evolution equation and the initial condition: 



2m 



(3) 



\im.KD[rf,T]r,,,Q\^5[rf~ri) . (4) 

tXO 

If a propagator containing two indistinguishable particles j and k has to be derived from a projection on the 
symmetric (anti-symmetric) irreducible representation of the permutation group has to be made for bosons (fermions): 

Ki[rf,T;r,,0]^^J2^'^^D[Pff,T;r,,0] , (5) 
p 

where P represents the permutations of the particles j,k and ^ — ±1 is chosen according to the particle statistics 
under consideration. The evolution equation for the projected propagator is given by the same eq. (^, because V| 
as well as the potential V (r) are invariant under all the permutations P. The initial condition obtained from the 
projection by taking the limit t J, is not the initial condition that is usually studied if one wants to solve a 
partial differential equation using probabilistic methods. [ p| In order to obtain an initial condition which is given by 

lim/f/[r/,T;r,,0] = 5(r/ -n) with r/,ri G Z)^ ® M^^^-^) ^ (g) 

rlO 

one has to restrict the configuration space M.^^ to a domain Z?! (g) R^(^~^). In this domain the identical-particle 
coordinates rj,rk are linearly ordered. In D2, for example, the ^-coordinates are ordered as follows: xj > Xk- 
Introducing a domain implies that one has to specify the boundary conditions; the specification appropriate for 
bosons or fermions moving freely or interacting harmonically has been proposed in Q and analyzed in 1 1[ . 
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1. Boundary conditions 



The ordering particles in a three-dimensional space can be performed by ordering the particles with respect to each 
of the three independent Cartesian directions. This option has been chosen in the previous investigations and led j^-^] 
to the construction of an ordered state-space D. The symmetric irreducible representation implies four combinations 
I = 0; 1,2,3. By definition, the {I = 0)-boundary condition causes reflection at the boundary dD in every direction 
X, y, z, whereas for Z = 1, 2, 3 the boundary dD reflects in x, y, z-direction and absorbs in the two other two Cartesian 
directions. For fermions, I — means absorption at dD in all three directions, whereas for I ~ 1 absorption occurs 
in the a;-direction and reflection in the y- and z-direction. Similarly, I — 2 and I — 3 correspond to the appropriate 
absorption and reflection conditions with respect to the y- and z-direction. It is clear that in order to study the 
evolution of the system in Euclidean time one has to give not only the starting point in the domain Z?! but also the 
boundary conditions I. The propagator has to contain the information that a state, initially characterized by Ti and 
I, will end after a Euclidean time lapse r in r^^ with boundary conditions I. The evolution in Euclidean time of this 
propagator can be studied from its extension to R'^^, which gives 

Ki[rfJ,T;rJ,0]^(^) J^^f'' Y.^^" T.^^' ^D[P.PvP^rf,T;rJ,0] , (7) 

^ Py Pz 

where Pj.,Py,Pz symbolize the permutations of the x, y, z-coordinates of the two indistinguishable particles j and k 
and where we imposed the initial condition 

\imKi[rf,l,T;n,i,0] = 5iiSirf-r,) , r/,r, e i?^ ® R^'^"') . (8) 
The boundary condition I determines the parity ^^^^ , , of the permutations Px,Py,Pz- 



2. Evolution equations 

A potential V (r) whose decomposition contains only one-dimensional irreducible representations with respect to 
permutations Qx,Qy,Qz of the coordinates is straightforward to analyze. Expanding the potential energy in the 
possible symmetry combinations, described above, one may write: 

3 3 
/'=0 ^ Qx Qy 

Consider now the following decomposition 

3 

2! 



X 



I) EEE^^^f^^r^ ViPxPyPzTf) Ku[P.PyP.rf,T-rJM 

Px Py P. 



/' PxQx PyQy P^Qz 

Introducing the permutations Si = QiPi with i — x,y, z, the quantity x can be written as: 
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X 



-Px S X Py S y Pz S z ^ ' 



—1 SP^^ Q P 

The permutation Pi and its inverse P^ have the same parity and therefore ' — We furthermore define 

= ipip . This leads immediately to: 

x = E^i(^/)^^[^/'^'^;^-^"'0]' (10) 
I 

where we introduced the notation: 
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Using eq. (||) for the propagator (Q) and the resuh obtained in ([l0|), the evolution equation becomes 

— Ki[rf, l,T;r,,i,0] = ^ V| Ki[rf,l,T;r,, 1, 0] - ;^ E ^^[^/' ^3 ' (1^) 



I 

where use has been made of the invariance of with respect to permutations of the coordinates x, y, z of the particle 
positions. 

The study of the non-diagonal interaction contribution (^) is the main purpose of the present section. In some 
cases, these non-diagonal terms are zero. For instance, if V (f ) can be written as Vg (r) = V (x) + V{y) + ¥{1), the 
invariance under the permutations of the particles implies invariance under permutations of the coordinates, which in 
turn then implies that Vii>{r) = 5ii'V{r). For harmonic interactions along the principal axes the many-body potential 
can be written as such a sum. The evolution equation is then diagonal in the quantum numbers Z, meaning that 
during the evolution the boundary conditions do not change. The computational feasibility and efficiency of this type 
of evolution equation has been demonstrated in jl] . 

3. Decomposition of the Potential 

The potential V{rj; fk) = Vi^i, . . . , r^, . . . , r^, . . . , rjv) is invariant under the permutation of any two particles, and 
in particular for interchanging fj with fk 

V{f,;fk)^V{fk;fj) ■ (13) 
For motion in 21?, one can decompose the potential as: 

V{xj,yj;xk,yk) ^ Vbh{xj,yj;xk,yk) + Vs{xj,yj;xk,yk) , (14) 

with 

Vbh{x j,yj]Xk,yk) = ^ [V{xj,yj;xk,yk) + V{xk,yj;Xj,yk) + V{xk,yk;xj,yj) + V{xj,yk;xk,yj)] 

and 

Vff{xj, yj;xk,yk) = ^ [V{xj,yj;xk,yk) - V{xk,yj; xj ,yk) + V{xk,yk; Xj ,yj) - V{xj,yk;xk,yj)] . 

Introducing the operator ax,(7y for interchanging the x,y -coordinates, and the operator ex,ey for leaving the x,y- 
coordinates unchanged, the decomposition of the potential can be written as: 

Vhh{xj, yj;xk,yk) = ^ (e^ + ctx) (e^ -I- ay) V{xj,yj;xk,yk) ■ 
For motion in 31? the same principle can be used leading to the decomposition: 

V{xj,yj,Zj]Xk,yk,Zk) = Vt,hh{xj,yj, Zj;xk,yk,zk) + V\,s{xj,y.j,Zj\Xk,yk,Zk) (15) 

+Vfhi{xj,yj,Zj;xk,yk, Zk) + Vsb{xj,yj, Zj;xk,yk, Zk) , 

where typically Vfhi[xj,yj,Zj;xk,yk, Zk) is given by: 

Vfbi{xj,yj, Zj-Xk,yk, Zk) = (e^ - a^) {cy + ay) (e^ - a^) V{xj,yj,Zj;xk,yk,Zk) (16) 

o 

and analogous definitions hold for the other matrix elements. It is interesting to note the similarity between this 
projection and quaternions: 



(ea; + a^) (cy + ay) (e^ + cr^) + {e^ + a^) (cy ~ ay) {e^ - a^) 
+ {<^x - a,j) {ey + ay) (cz - a^) + (e^ - a^) {cy - ay) (e^ -f- a^) 



^xyz — g 

Using this projection the formal decomposition (InJ) is achieved explicitly. 



(17) 
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B. The Feynman-Kac functional 



Given that the propagator indeed satisfies the evolution equation (12), it remains to be shown that a process can 
be found to generate the sample paths and a Feynman-Kac functional. 



1. Positivity and the process 



To start the construction of the generator of the process consider the following identity: 

I Sy 



21 

V{r) 



This is a direct consequence of eq. (^. Assume further, that the domain is chosen in such a way that V^j (r) < 
for / 7^ r . This allows to rewrite eq. (f2h as follows: 



— Ki[rf,l,T-rJ,Q\ ^—V%Ki[rf,l,T-rJM ~ -V{rf) Ki[rf,l,T-rJ,Q] (19) 

-Wi (Tf) Ki[rfJ,T;f,,i,0] + ^ M^^j (Tf) Ki[rf, l,T;r,, 1,0] , 

I 

where Wjj (Tf) denotes the non-diagonal part —jVii{rf)iil^l: 

1 

li 



Wu{rf) = -\va{rf){l-5,-;) , (20) 



and the quantity Wi{rf) stands for: 

W^;(^/) = E^»(^/) (21) 

In this situation, the Wjj (rj) can be interpreted as the rate that a walker in the position r f will change its boundary 
condition from liol. This means that the following terms of eq. (|9|) generate a diffusion together with a Markov 
process on the four possible boundary conditions: 

|-7^0[r/,Z,r;r„f,0] = A ii:0[r^, Z,r; r,, f , 0] -f ^ l^^j (r^) if,"[r/, r,T; r,, f , 0] (22) 

I 

-Wi{rf)K''i[rfA,T-rJ,% 

where K'j\ff, l,T;ri, 1 , 0] satisfies a backward Kolmogorov equation, and there for it is the transition probability of a 
compound stochastic process pO| that generates the sample paths in our state space. 



2. Positivity and the functional 

Above, we assumed that the domain ZJf can be chosen such that VJy (r) < for Z ^ I. This assumption is 
necessary, because otherwise the introduced VF-matrix would not represent a continuous-time Markov process. From 
its construction it is clear that it depends on the potential. Therefore a general analysis for any potential is beyond the 
scope of the present communication. Assuming that the system under consideration has a potential with this property 
(or can be transformed by a unitary transformation to a system with such a potential), the evolution ( p2|) describes 
a diffusion triggering a continuous-time Markov chain that swithces boundary conditions in a Wiener-Poisson space. 
The switches do not induce discontinuities in the sample paths but allow for a change in boundary condition according 
to the rates defined by the H^-matrix. The sign condition on the projections of the potential is necessary for the 
interpretation of the non diagonal elements as rates. The evolution equation based on this assumption describes a 
stochastic process in the strict mathematical sense. It can therefore be simulated sign-problem-free. The propagator 
Ki[rf,l,T;ri,i,0] with the initial condition (||) can be written as a Feynman-Kac functional 
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where use has been made of the identity ( ]l8| ) to make exphcit that the potential is diagonal in the label I. 

The sample paths are generated by a 3A^-dimensional Brownian motion {i? (r) ;t > O} with — fiT/m and ab- 
sorbing and reflecting boundary conditions as defined by the index I. After the transition, the boundary conditions 
are adapted to agree with those from the Markov chain {Lt]t > 0} where the (normalized) quantities Vjj then deter- 
mine the boundary-condition transition rates I ^ I . Repeated application of this procedure for all the subprocesses 
starting with boundary conditions I generates the sample paths of the diffusion process of — 2 distinguishable and 
two indistinguishable particles in three dimensions. The stochastic integration over the Feynman-Kac functional is 
then done by killing and branching. 



=i\ exp 



(23) 



III. DIFFUSION ON A REDUCED STATE SPACE D 



The positivity of the rates Wfi{rf) is a necessary ingredient to guarantee the interpretation of the Euclidean time 
Schrodinger equation as a compound diffusion equation with killing and branching in a Wiener-Poisson space. For 
some potentials, the condition of negative definite V'jj(r/) may not be satisfied on the state space 

cases, the definition of an alternative (reduced) state space is required in order to find purely positive transition rates 
Wij{ff). The symmetry of the potential effects the structure of the adapted state space, which therewith is specific 



for certain classes of applications, 
be considered, where 



To illustrate the construction principle, the reduced state space D 



^2 — 



Xj,xk,yj,yk,zj,zk) Xj > \xk\,yj > \yk\,Zj > \zk 



1} 



will 



(24) 



This choice will be of relevance for the probabilistic sampling of the real ground-state wave function of ortho-helium. 
The procedures for the construction of an ordered state space Xj > Xk, yj > yk, Zj > Zk, reported in |^-|^, alone do 
not induce the ordering on ID2, and an additional set of boundary conditions must be supplied to impose the desired 
ordering. In that framework, the separation of the multi-dimensional diffusion process into one-dimensional diffusion 
processes will again prove useful. 



A. Free Diffusion of Indistinguishable Particles on 

The free diffusion process of two indistinguishable particles on a line has been introduced in as the diffusion 
process of two distinguishable particles on the ordered state space with absorption (for fcrmions j or reflection (for 
bosons) at the boundary ■ The construction of the free diffusion process on the reduced state space requires the 
extension of the stochastic process. (Anti)symmetrization again plays the major role in that approach: the diffusion 
process of free indistinguishable particles on the state space can be mapped onto D2 by the construction of 
reflecting and absorbing boundaries. The starting point for the derivation of the required processes is the separation 
of the state space into two non-intersecting subspaces D2 and PxD\ of equal volume: D\ = D\\J PxD\. The 
mapping between both subspaces is realized by the permutation operator which models reflection at xi = —X2, 
i.e.: Px{xi, X2)^ — {~X2, —xi)'^ . Then, the transition probability density p{^h{x,T;x') that the system of two free 
fermions/bosons evolves from x' to x during the Euclidean time interval r, can be written as 

Pi^h{x,T;x') ^ ^{pI]^(x,t;x') + pI^{x,t;x')) , (25) 

where we defined 

Pf b(x, t;x') = pthix, t;x') - p^hiPxX, t;x') , p^b(^i t;x') = pf,b(x, t;x') + pib{PxX,T;x') . 

Each of the transition probability densities p{^hix,T;x') consists of two density matrices, one — p^y^{x,T;x') — for 
which the boundary xi = —X2 acts reflecting, and one — pf^{x,T;x') — which induces absorption at xi = —X2- With 
an appropriate initial choice x' G D^, free fermion/boson diffusion to any point x E D\ — and thus also to any 
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X G — can be projected to the reduced state space D\ by means of reflecting and absorbing state-space boundaries. 
Each of the elementary density matrices p^^{x,T;x') and p\^{XtT]x') constitutes a diffusion process with the initial 
condition 



lim /9f b(a^j T^x ) — lim Pf b(a;, t\x ) — 6{x — x ) . (26) 

tIO ' tJ.0 ' 

To proof this statement, pf{x,T]x') and Pf{x,T;x') must be shown to satisfy probability conservation and the semi- 
group property. Both properties are derived taking into account the Markov properties for the free fermion/boson 
diffusion density pf i^{x,T;x') defined on the state space D2 with absorbing/reflecting boundaries. Probability con- 
servation for the two-boson density p^{x, t;x') follows from 

/ dx p^,{x,t;x') y^,Ph{PxX,T;x') = | / dx pI{x,t;x') ^ [ {^Ph{x,T;x') + pb{PxX,T;x')^ = 1. 

As argued in j^, the free fermion diffusion density p{{x,T;x') is not automatically normalized on the state space Dj- 
In contrast, normalization of the free fermion diffusion density p{{x,T]x') requires the introduction of the absorbing 
boundary state j|] which comprises all transitions to the outside of the state space . Also in the present analysis, the 
construction of a boundary state [ pl| ensures the normalization of the probability densities pI^(x,t;x'), p^{x,t;x') 
and Pf{x, t;x'). The density pl,{x, t;x'), for instance, associates an absorbing boundary xi — —X2 to bear in mind the 
probability for distinguishable-particle propagation to a;i < —X2- The semigroup properties of Pf(;E,T;x ), p^{x,t;x ), 
p{,{x,t;x ) and p^{x,T;x ) are derived by unfolding to the state space D^. Using their symmetry properties, they 
can be retraced to the semigroup properties of pf{x,T;x ) and pb{x,T]x ). E.g. 

/ dx p\{x' ,T]x) p\{x,(t;x") ^\ I dx^^ p\{x' ,t;Pxx) p\{Pj:X,a;x") 

= ^ l^x (^pf{x' ,t;x) - pi{x' ,t; Pj:x)^ (^p f{x,a;x") - pi{PxX,a;x")'j 

= / dx (^p{{x',T]x) p{{x,a]x'') - pf{x',T]x) p{{x,a-,Pj:x'')j 
= /9f(x', r + cr; x") ~ pf(x' ,t + a; P^x") = p\{x' ,t + a;x") . 

As a consequence, the free diffusion process X{^b{T) of two indistinguishable particles on a line is characterized by the 
combination of two adapted Wiener processes defined on the reduced state space : 



Xhir) = i (^^^(t) +xI{t)) for bosons, and 



Xf{T) = i (^Xf{T) + X|(t)) for fermions 



Since p\!^f{x,T;x ) and p\^f{x,T;x ) are orthogonal, e.g. 

r r 

dx pI{x\t]x) p^(x,t;x'') ~ / dx pI{x\t;x) p^{x,T-,x'') + dx Pf{x' ,t\x) p^{x,T-,x'') 

D\ Jd}, Jp^Di 



/ dx [pI{x',t;x) p\{x,t;x") - Pf{x',T;x) p}{x,t;x")] =0, 



each of the densities p^i(x,T]x ), pl^f{x,T;x ) evolves independently and generates the corresponding Brownian 
motion. We have used in the upper and the lower indices the same letter to characterize the boundary condition, this 
is customarily done in supersymmetric models it should be noted that these letters for b in the upper indices 
indicate wether there is absorption or reflection on the boundary introduced by the (broken) symmetry requirements 
for the potential and they are independent of the f or b in the under indices introduced by the statistics of the particles. 
In the sampling procedures, each of the processes X^'^^t) is constructed from the diffusion process -'^d(''') of two 

distinguishable particles taking into account the corresponding boundary conditions: the process xI{t), for instance, 
means reflection at the hypersurface xi — X2 and absorption at xi — —X2- For each of the composite processes, 
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absorption at the boundary defines a characteristic Erst exit time t^jji Q which specifies the adapted stochastic 
process, e.g. 

Xdir) for T < Tg^i 
^dCrgfti) for T>Tg5i 



In fig. 1, the feasibility of the developed sampling technique is indicated by graphical comparison of the numerically 
sampled and rigorous transition probability densities Pf for a given initial x' denoted by the crosses. 



B. Free Fermion and Free Boson Diffusion on Z)? 



We now formulate the diffusion process of two free indistinguishable particles in three dimensions on the reduced 
state space by extending the symmetry analysis performed above. Following the basic principles outlined in 
the free diffusion of indistinguishable particles can be traced back to one-dimensional free diffusion processes. In 
that framework, the transition probabilities Xf'^{T) represent the elementary subprocesses combined to constitute the 
overall free diffusion process in three dimensions. For free identical fermions in three dimensions, the density matrix 
can be expressed as 



p(r,r;r') = ^ 



[pI{x,t;x') +p^(x,t;x')] [pi{y,T;y') + pl{y,T;y')] U (z,t;z') + pI(z,t;z')] 



[pi 



Pt 



Pf\ [pi +pi] + [pi +pi] [pi 

whereas for the corresponding boson system the density matrix evaluates to 



[pf + Pf] [Pi+ Pf] [pI+pY. 



p{r,T;r') = ^ 



[piix,' 
+ [Pi^ 



x') + pI(x,t;x')] [pi{y,T;y') + pl{y,T;y')][pi{z,T;z') + pI{z,t;z')] 



■ Pn [Pi + Pr\ [pi + Ph\ + LPt + Ph\ [Pi + Pf. 
In this representation, the density matrix for identical particles contains 32 terms. Since each term factorizes into a 
product of three one-dimensional orthogonal density matrices, each of the 32 elementary three-dimensional density 
matrices satisfies the Markov properties and generates a Wiener subprocess on the reduced state space Z?2- Fo^' 
instance, the subprocess X^I{t) 



[P'r 



P\ 



[P\ 



■P\] 



[pi 



pI 



[p\ 



p)] 



^fbb 



{T)^Xl{T)®Yl^{T)®Zi{T) , 



can be sampled as the diffusion process of two free distinguishable particles restricted to the reduced state space by 
reflection at the boundaries yi = y2 and zi = Z2 and absorption at Xi — X2, Xi — —X2, yi — — 2/2 and zi = —Z2- The 
combination of all the 64 Wiener subprocesses gives rise to the free diffusion process of two indistinguishable particles 
in three dimensions defined on Z)| . The symmetry properties of the elementary processes allow furthermore to unfold 
the resulting probability density into the configuration space R^. 



C. Evolution Equations and Decomposition of the Potential 

The derivation of the evolution equation for the propagator proceeds along the same lines as for the two-dimensional 
case, as discussed in Sec. II. A. 2. Without going into detail, we only mention the result 

— Ki[rfJ,l',T-rJj',Q] ^—V^Ki[ff,l,l',T;nJ,i',0] -jJZ.^a'^^s) Ki[rf,l,l\T;n,iJ',0], 

i~i' 

where 

Ki[rfjx,T-rj:i'M = (^)'e^/'^ E^^ E?;'^ E^r^ E^f" E^r^^i^i^-^?^^-^-^?^^^^/'^;^-^"'^] ■ 

^ p. Py p. Py P^ 

The components of the potential are given by 




it should be noted that if the range of the label I used to denoted the statistics is extented in such a way that it also 
indicates the boundary conditions on the reduced domain then the evolution equation for the propagator derived here 
takes the same form as the equation (|T^). 
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IV. AN ALGORITHMIC APPROACH ON MODELS 



In the previous section, we introduced a novel representation of the Euchdean-timc propagator of systems containing 
two indistinguishable particles. This representation allows to sample the propagator sign-problem- free. We now 
develop a numerically feasible sign-problem-free implementation scheme. Its efficiency in numerical practice is analyzed 
in applications to model systems. In all the cases, we use atomic units h — m — e — 1. Energy and Euclidean time 
are measured in Hartrees (H) and H~^. 



A. General Structure 



The key elements of the many-body diffusion algorithm (MBDA) for the symmetry problems considered here are 
(i) an appropriately adapted free diffusion step, (ii) a routine to model subprocess transitions, and (iii) a branching 
and killing procedure. Whilst some elements of the algorithm resemble Diffusion Monte Carlo, we stress that the 
underlying principle of interdependent subprocess evolution is conceptually new. We are also aware that the efficiency 
of the presented algorithm can be drastically enhanced by importance sampling methods. Although we can show that 
this expectation holds true, we do not dwell on these refinements, for major emphasis is on the realization of symmetry 
arguments. In the limit of infinitely long evolution r — > oo, the stepwise application of the Euclidean time propagator 
on a given sample leads in principle to the system's lowest eigenstate compatible with the symmetry restrictions 
made. In numerical practice, this evolution occurs in discrete Euclidean time steps At. Since the assumption of such 
time steps is only exact for At —^ 0, the algorithm in principle suffers from a systematic time-step error. The latter, 
however, could be controlled in the considered applications. The main advantages of the developed approach are the 
simplicity and convergence properties. In particular, it allows for the transparent implementation of reflecting and 
absorbing boundaries to model the subprocesses derived above. Let us discuss the algorithm in some more detail. 
Having generated an initial population of walkers located at positions in the appropriate state space, each valid 
walker makes a free distinguishable particle move to {ff)i- This is achieved by the construction of Gaussian deviates 

with variance V At and mean f/. If a walker hits the state-space boundary by free distinguishable particle diffusion 
it is either reflected or absorbed depending on its characteristic boundary conditions [|l|. In numerical practice, the 
use of non-zero Euclidean time steps At prevents the proper construction of reflecting and absorbing boundaries. 
In the mathematical literature, free diffusion in the presence of an absorbing boundary is realized by the so-called 
Skohorod construction The latter reflects walker trajectories without changing the "momentum" of the walker. 
This is in contrast to the reflection experienced by a classical object at a hard wall. In the numerical procedures, the 
final reflected position r'j is obtained by reflection of the corresponding f'j. The systematic error for inappropriate 
reflection vanishes for At 0; for the utilized time intervals At no significant error were observed. The efficient 
implementation of fermionlike diffusion, in contrast, calls for numerical crossing-recrossing corrections. For two 
identical fermions j, fc on a line the corrections for an absorbing boundary Xj = Xk have been discussed in Q . Based 
on the principle of images, we showed that an apparently valid walker move from (x'j,x'j^) to {xj,Xk), both points inside 
the state space, has to be rejected with the probability exp(— (xj — Xk){x'j — x'i^)/At). The symmetry properties of the 
potentials discussed in this work require the generation of interdependent subprocesses. In succession to the adapted 
free diffusion step, a walker associated with the boundary condition of the subprocess I may change to the boundary 
condition of a subprocess I' with the probability pi^i' {ff). Its actual structure is determined by the physical situation, 
e.g. the symmetry properties of both the potential and the eigenstate under consideration. The local dependence of 
the probabilities pi^ii{ff) governs the relevance of subprocess transitions at a given position ff. The structure of 
pi-^ii{ff) is fairly flexible as long as it harmonizes with the killing and branching procedure. For instance, pi^ii{ff) 
could be decoupled from the different kinds of boundary conditions. Then each subprocess transition occurs with 
equal probability , where ni denotes the total number of subprocesses. Accordingly, both the normalization rates 
and the probability for branching and killing become dependent on the subprocess transition performed on the walker. 
In what follows, we assume subprocess transition probabilities adapted to the branching and killing procedure due 
to the potential. The third step, branching and killing, realizes the presence of the potential by the implementation 

/ Ar _ ^ ^(^^)=^ 

of the exponential path-integral weights w{r',r) — exp — / (^('?)) 







, where i?(<?) denotes the paths 

R(0}=f' 

generated by the adapted free diffusion process. The principal lack of continuous trajectories in discrete time-step 
evolution procedures induces the necessity for feasible approximations to w{f', f). For sufficiently smooth potentials, 
the Suzuki- Trotter formula offers a reliable alternative. Then only the values of the potential at the initial and final 
points of a move are taken into account: w{f',f) = exp{—{V(f') + V{f))AT/2}. Especially for singular potentials. 



improved approximations (for a detailed discussion see e.g. |7| and references therein) are recommended. Here, we 



ai 
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choose the semiclassical approximation suggested in w(f f) = exp{— d£, V{f' + ^{f — r'))}. A branching and 
killing technique serves to efficiently adapt the walkers to the weights w(f',f) by an acceptance/rejection procedure. 
To conserve the size of the walker population to within statistical fluctuations, an appropriate reference energy iS^ef 
has been supplied to yield improved weig hts w(r',f)e^-f'^^. In H, we derived an estimator for the ground-state 
energy Eq. The key elements of this estimator are the potential average {V) over the state space D and a surface 
term (j)^^, triggered by the non-zero gradient of the fermionlike subdensities at absorbing boundaries. The estimate 
is easily generalized to an ensemble of interdependent subdensities 5'(r, 1,t) 



Ef) — lim E. 



V = lim y 



E 



'V) 



Alternatively, the growth estimates 

Eq — lim Et- ~ — lim — In / drr [ dTi'S^ K[rf,l,T;ri,0] , 

T— *oo r^oo T J J — ^ 

may serve to predict the ground-state energy. Based on this general algorithmic structure, the ground-state wave 
function of several model systems is investigated below. The different types of potentials make a specific analysis of 
the appropriate state space necessary. 



B. Double well and Mexican Hat 



To start with, the symmetry principles underlying the many-body diffusion algorithm are illustrated for a quantum- 
mechanical particle in a double-well or a Mexican Hat potential. For these models, the utilization of (anti)symmetric 
subprocesses is shown to provide a suitable technique for the numerical generation of the real ground-state wave 
function. Our study aims at revealing the role of symmetry arguments to sample a probability density on an ap- 
propriate part of the configuration space. Feasibility of the symmetry concept has been tested independently with 
standard Monte Carlo techniques. Consider a single particle in a one-dimensional double-well potential, described by 
the Hamiltonian 

H=-^-^ + V{x) with V{x)=^-^^ - + bx\ (27) 

where a and b represent parameters. Fig. 2 depicts the potential values for two different parameter choices: a = 2. = 
and a = 2,b = 0.25 in the following referred to as symmetric and asymmetric potential. According to the rules 
outlined above, the Euclidean time propagator K{x,t;x') associated with the Hamiltonian ( p7| ) can be formulated 
as a sum of four subpropagators K{x,t, j; x' , j) with j,j € {0, 1} denoting symmetric resp. antisymmetric behavior 
under the transformation : x ^ —x. This means that K{x, j,t; x' ,0) and K{x,j,t;x', 1) are initially symmetric 
and antisymmetric with respect to the boundary x = which separates the two half-spaces D~ = {x\x < 0} and 
= {x\x > 0}. The (anti)symmetrized potentials V^(a:) and Va{x) are derived as 

rr, ^ TW N V{X)+V{~X) x'^{x'^~2a) , , , . V{x)-V{-X) ,3 , 11 / / 



The matrix W (cfr. (|20[)) is calculated as W,,j i = — (1 — Sj^/) with the total rate represented by Wj = — Vs.. The 
Euclidean-time Schrodinger equation [cfr. (119)] can thus be written 



2 



1 



d 1 

— K[x, j,t; x',j, 0]^-—K[x, j,t; x', j , 0] - Vsix) Ki [x, j,t; x',j, 0] - }_J V,{x) K[x, / ,r; x' , j , 0] (1 - <5,,y ) . (28) 

^ f=o 

Eq. ( |2^ ) allows for a change of the boundary conditions of the propagator K[x,j,T;x',j,0] during its evolution in 
Euclidean time. Written as a Feynman-Kac functional, the propagator K[x, j,t; x' , j,0] can be interpreted as a sum 
of two interdependent distinguishable-particle diffusion subprocesses. The numerical realization of the subprocesses 
J occurs by means of ensembles (a;:^)i=i.„^ of n-j walkers. Each of the rij walkers satisfies the boundary conditions 
for the subprocess j and during the evolution a walker may change its subprocess j to j according to the normalized 
transition amplitudes Pj^f{ff). To formulate a probability-based evolution scheme it is therefore required that 
Pj-,y{ff) remains positive. Since 
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^exp i-rVfjirf)) 
3=0 



2 exp 



{x^ - af 



[cosh (- T 6 x^) Sf^j + sinh (- rbx^) (1 - Sy^j)] 



the half-space D~ represents the appropriate state space for probabilistic sampling. Walker transitions to x > are 
prevented by the construction of a reflecting (for j=0) or an absorbing (for j=l) boundary. The major algorithmic 
steps for the numerical evolution procedure are enlisted as follows: 

1. Generate two distinct walker sets (a;^)i=i,„^ , j=0, 1. The initial non-zero walker populations rij are chosen 
randomly. 

2. Calculate for each walker the subprocess-transition probabilities Pa,{x) = [l — exp (— 2AT6a;'^)] /2 and apply 
them by comparison with uniform pseudo-random numbers r/ e [0,1]: if r? < Pf(x) change the boundary 
condition associated to the walker, else leave the condition unchanged. 

3. Make an evolution step, as indicated in Sec. A. 

4. Return to step 1. until the density has equilibrated; after a certain amount of repetitions, calculate the potential 
average and the walker outflow to get an estimate for the ground-state energy. 

5. To obtain an image of the sampled ground-state wave function, record the walker positions in small discrete 
bins. 

The results obtained for the ground-state energy Eq of the double- well system arc shown in table 1. To test the 
algorithm, the estimates Eq^^^-^ have been compared with the numerical outcome of a program with evolution 

on the full configuration space. For both the symmetric and the asymmetric potential, the results agree within 
the estimated standard deviation a. It is illuminating to examine the population sizes of the two subproccsscs j. 
For the symmetric potential the probability Pa.{x) is zero and independent subprocesses are simulated. Due to the 
higher eigenenergy related to the fermion subprocess X[x, 0,t; x', j, 0], the antisymmetric process rapidly fades out 
and the ground state wave function is completely described by the symmetric process. For the asymmetric potential 
in contrast, the antisymmetric subprocess is of crucial importance, since the asymmetric ground-state wave function 
cannot be simply generated by a symmetric process. The faster decay of the antisymmetric density is counteracted 
by a net walker transition from the symmetric to the antisymmetric population. In equilibrium, almost half of the 
population belongs to the antisymmetric density. 
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6 = 
b = 0.25 


-0.2996 
-1.0251 


0.0006 
0.0012 


-0.2999 
-1.0251 


0.0006 
0.0011 


100 % 
53.7 % 


% 
46.3 % 



Table 1 



This example for a particle in one dimension can be transparently generalized to arbitrary spatial dimensions. 
Consider, e.g., a particle in two dimensions exposed to the Mexican-Hat potential 



1 T.. ^ 



with V{x,y) = 



(a;2 - af 



+ bx^ 



(29) 



The two parameter sets a = 2, 6 = and a = 2,b = 0.25 were chosen as in one dimension to represent a symmetric 



and an asymmetric potential, as depicted in fig. 2. Consider the two operators x ^ —x and Py : y 



-y and 



introduce the two boundaries a; = and y = 0. The construction of reflecting and absorbing boundary conditions 
at these two boundaries calls for the deflnition of 16 Euclidean-time propagators K{x,T,j;x',j) and four elementary 
states ^ss, *sa, ^'as and ^aa with their indices characterizing their properties under transformation with P^ and Py, 
respectively. By imposing reflecting and absorbing boundary conditions the configuration space is split into four 
subspaces with positive or negative x- or y-coordinates. Although the system's evolution can be folded to any of 
these four subspaces, the analysis of the subprocess-transition amplitudes for the Mexican-Hat potential reveals that 
only two subspaces are associated with positive definite sampling. A sign analysis of the four subprocess-transition 
amplitudes Pss{rf), Panir/), Pas{rf) and Pai,{ff) 



Pas{x) 



I - exp {-2ATbx^) 



Pssix) 



1 + exp (-2At6x3) 



Psnix) = , Pa^{x) = 



indicates that they all may be interpreted as probabilities if the evolution is confined to either the subspace Dt = 
{x < 0,y > 0} or DZ = {x < 0,y < 0}. The fact that Psa{x) = and Psia.{x) = induces two independent sets of 
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subprocesses, namely [K{x,T,ss;x',j), K{x,T,a.s; x' , j)] and [K{x,T,sa;x' K{x,T,aia,;x',j)]. Only those subprocess 
transitions are allowed which conserve the parity under Py. The formulation of the respective MBDA yields results 
for the ground-state energy which are in good agreement with the results achieved without the introduction of 
boundaries (see table 2). Since the characteristics of the model ( p9|) allow for a distinguishable-particle treatment, 
the boundary-free algorithm reliably simulates the properties under investigation. The comparison of the estimated 
ground-state energies demonstrates the efficiency of the MBDA. Both energies and E^^^^ have been obtained 

with comparable numerical effort and standard deviations. The last four columns of table 2 show the relevance of the 
four states ^'ss , ^'sa, ^as and ^'aa with as measure of their relative population sizes. Since Psa,{x) = and Paii{x) = 0, 
the states Vl/ga and \l/aa fade out. In addition, for the symmetric potential Pas{x) = 0. This means that also ^'as is 
irrelevant for the sampling of the ground-state wave function. The remaining state is symmetric under both P^; 
and Py. Clearly, evolution in the presence of the asymmetric potential cannot be simulated with a symmetric state 
only; the algorithm ensures the occurrence of the "ifas state by a steady transformation of walkers associated with the 
subprocess K{x,t,ss; x' , j) to walkers propagating according to K{x,T,e.s;x',j). 
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E*as/E* 


E*sa/E* 


E*aa/E* 


6 = 
6 = 0.25 


-0.19865 
-0.56129 


0.00042 
0.00054 


-0.19860 
-0.56115 


0.00037 
0.00050 


100 % 
60.5 % 


% 
39.5 % 


% 
% 


% 
% 



C. Anisotropic Harmonic Fermions 



A system of two non-interacting anisotropic identical fcrmion oscillators in three dimensions serves as a transparent 
testing ground for indistinguishable particles. We considered a system with the Hamiltonian 



H-^ + ^{xl+xl) + '-f{yl + yl) + '-f{zl + zl) , 



(30) 



with ground-state energy 

We now consider a rotated reference frame as the result of from the rotation by the two Euler angles 4> and 9 |l^] 



R = 



cos (p 
— sin (j) cos ( 
sin 6 sin 6 



sin (j) 
cos cos 9 
— cos 6 sin 9 cos ( 



The objectives are clear: the symmetry of the potential under the interchange of Cartesian particle coordinates is 
avoided without losing the advantage of comparison with a rigorous analytical solution. A proper implementation of 
the propagator with the MBDA requires positive definite subprocess transition amplitudes 



Es. Es„ Es. ef^CrCf^- exp i-rViS^SyS. Vf)) 



E7E5. E. E5. "Cf^ exp {^rViS^SyS. Vf)) ^^^^C 



(32) 



where we used the notation introduced in Sec. II. If the numerator of cq. ( |32D is positive definite, Pijirf) may 
be interpreted as the probability for a subprocess transition I —t I. Since the sign oi Pij{ff) is determined by the 
parity of the propagator if [r/,Z,AT; r^, Z,0] and the boundary condition I under Cartesian-coordinate permutation, 
four different subprocess-transition amplitudes occur. Their classification according to their symmetry properties 
suggests the following notation: Phhhijf), PSh{rf), P[hi(ff) and pbff('^/)- The probability psh(ff), e.g., stands for 
transitions between subprocesses with even parity in the z-coordinate and odd parity in the x, y-coordinates, for 
instance i4r[r/, fbb,Ar;ri,bfb,0]. While the amplitude phhh{r) is always positive, positivity of the three remaining 
amplitudes is imposed (in first order of the Euclidean time step At) by the following conditions on the domain and 
the parameters fix, ^y, and 9: 



sin (j) cos (j) - 
cos cj) sin 9 cos 9 
— sin 6 sin 9 cos ( 



nl cos^ ( 



^ rtlicos^ 9 - 1)) {xi 
{yi- y2) {zi ~ Z2) 

fll) {xi - X2) {zi - Z2) 



X2) iyi - y2) 




(33) 
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For the present illustration, the evolution has been restricted to the state space -D| with a particular set of parameters: 

xi>X2,yi>y2,zi>z2,n,^4., ily ^ 3 , ^ 2 , <j> ^ 130° , = 50° . (34) 

Employing samples of approximately 20,000 walkers and Euclidean time steps of 0.001 H^^, the MBDA leads to the 
numerical results visualized in fig. 5. The ground-state energy estimates (see left-hand side of fig. 5) oscillate in 
a small interval around their average of 10.9987±0.0056 H. They estimate the rigorous ground-state energy of 11.0 
H with a statistical inaccuracy of about half a per mil. After numerical thermalization, the population of each of 
the four subprocesses {f,b,b}, {b,f,b}, {b,b,f} and {f,f,f} remains approximately constant (see right-hand side of fig. 
5). In particular, contributions from the {f,f,f}-subprocess prove indispensable to correctly generate the propagator. 
Our numerical investigations manifested the occurrence of a net walker flow from the subdensities ^I'fbb) *bfb, ^bbf 
to '^gf. Without the inclusion of subprocess interdependencies, the subprocess K\r f ,Si,AT;ri,l,0] would fade out 
exponentially in evolution time. 



D. Application to Ortho Helium 

To illustrate the feasibility of the reduced state space concept, we also used the MBDA to calculate the ground- 
state energy of ortho-helium. In the triplet states of ortho-helium, the two electrons carry equal spin. The occurrence 
of additional repulsive forces due to exchange distinctly influences the energy spectrum, so that the simulation of 
the lowest ortho-helium state requires the accurate inclusion of Fermi statistics. In that context, the challenge of 
numerical implementation shows twofold. First, the lowest energy level should be reliably predicted and second, 
the converged energy estimates should be stable in time. In particular the second aspect concerns the fermion 
sign problem: the MBDA is demonstrated to exactly filter out possible interfering symmetric contributions. The 
general evolution scheme on the reduced state space D f relies on the construction of 32^ interdependent subprocesses 
K\ff,l,l',AT;ri, l, l',0]. Here, each of the 1,1 may be one of the four boundary conditions fff,fbb,bfb or bbf, whereas 
I', I' run over the eight conditions bbb,fbb,..,fff. Symmetry arguments allow to simplify the computational procedure 
by the prior rejection of certain classes of subprocesses. As the permutation operator P has the same effect on the 
ortho-helium Hamiltonian as the permutation operator P, the ground-state wave function consists of substates with 
equal parity under P and P. Consequently, the simulation of the ortho-helium ground-state wave function on the 
reduced state space D2 is accomplished by the construction of only 16 subprocesses with I = I, I' = I' specified as 
fff,fbb,bfb or bbf. Although both - the general and the simplified - procedure in principle lead to the exact densities, 
the restriction to four subprocesses considerably improves the convergence of the algorithm. The next step in the 
formulation of an evolution scheme for the ortho-helium ground state is the consideration of the subprocess-transition 
amplitudes 

pi'j ^ - ,^ ^ ^s.,Sy,s.,s^.,Sy.,sJi ^i' ^i' g ^1 y ^^^^ 

2^1,1'. S^,Sy,S,,S^,Sy,S,'il M ^/ ?i' ?i' e ^ " t.j t;j t;j, t-, 

controlling the interplay between the four single subprocesses. In order to derive the subspaces associated with 
positive sums of exponentials e'^^^^^^y^'^'^^y^'' '^f\ we focus on the short-time limits (1 — TV{SxSySzSxSySzr j)) 
of the exponentials, which contain all the essential symmetry characteristics necessary for the sign analysis of the 
amplitudes (^5|). Furthermore, it proves useful to separate the potential V{f) in two parts Vi(r) and V2(r) with 

^i(r) = ^ H- ^ , V2{r) = — ^ , V{r) = V^{r) + V2{r) . 
\ri\ \r2\ |ri-r2| 

To illustrate the need of a reduced state space, consider the corresponding many-body diffusion process on the 
state space Df. The amplitudes Phhh^j), Pffb(^), Pfhf^r) ^'^d Phs{f) control the transitions between the 16 fermion 
subprocesses if[r/, Z,At; r^, /,0]. A detailed investigation of the amplitudes, however, discloses that some of them may 
become negative on the state space 13 1. For instance, considering pshir), one derives 

Pffb(r) e^^(^) ^ {- f-{- fy{l-V,{P.xPyP.r)) 
4r 



dt 



13 



which is positive definite on a subdomain of the state space, namely for \xi\ > \x2\ and \yi\ >\y2\- Similar restrictions 
follow from pfbf(r) and pbff(^), whereas pbbb(^) is a sum of only positive elements. Strictly speaking, utilizing the 
state space as the domain for evolution, the MBDA would incorrectly simulate the ortho-helium propagator. It is 
at that point that the reduced state space l3f as defined by ( p^ comes into play. By appropriate restriction of the 
evolution to the reduced state space, the subprocess-transition amplitudes all remain positive and their probabilistic 
implementation is justified. In what follows, we take over the results from Sec. Ill and adapt the general scheme to 
the ortho-helium system. 

Continuing with the investigation of the amplitudes p(r), let us first focus on the contributions of Vi{r). In Vi{r), the 
coordinates xi, . . . , Z2 are entangled such that the transformations P and P have the same effects; for our purpose, 
no distinction has to be made between P and P. 

With regard to amplitudes with identical sub- and superindices per coordinate, one obtains 



P P P P P P 



The remaining 28 amplitudes have in common that for at least one coordinate the permutations P and P induce 
different parity. Each of these subprocess-transition amplitudes involves a sum composed of two parts of equal absolute 
value but opposite sign. Consequently, they are irrelevant for the evolution procedure. 

The potential V2{r) proves invariant under interchange of two coordinates. Its symmetry properties entail that only 
those amplitudes p{r) contribute which in each direction involve equal positive (bosonlike) parity under P and P, i.e. 



Pbbb(^)|y, « E cM-rP.PyP.P.PyP.r)^^exp{~rV2ir)) . (36) 



Obviously, Pbbb(f) is positive on the reduced state space D\. By recombination of both potentials Vi{r) and V2(r), 
the total evolution scheme is established. Its proper implementation guarantees accurate probabilistic sampling 
of the ortho- helium propagator on the reduced state space Df- The sampled density consists of 4 subdensities 
*fbb(^' t), ^''^b^(r, r), *bbf(f , r) and ^'|f(r, r). Due to the symmetry properties of the ortho-hehum potential V{r) = 
Vi (r) -I- V2 (r) , certain classes of subprocess transitions are forbidden; each of the subprocesses undergoes one out of 
four possible transitions ruled by the probabilities pg[^(r), p^j(r), Pbf(r) and Pbbb(^)- During a time step At, the 
subdensity ^'^[^(r, t), for instance, branches into the subdensities \E'b[j^(r, t -|- At), ^'[^{^[(r, t -|- At), 5'ff(r, t -|- At) 
and ^fbb(^, T + At) according to the probabilities p|{^(r), Pf|^f(r), Pbg('f) and Pbbb(^)- Notice that the analysis of the 
subprocess transition amplitudes is consistent with the restriction to four subdensities in the sense that no transitions 
occur to subdensities different from the four specified ones. Fig. 6 illustrates the energy estimates for the lowest ortho- 
helium state achieved with a pre-thermalized sample. The numerical average -2.1744±0.0012 H is in good agreement 
with the numerical energy estimate of -2.175229 H found in |lj]. Systematic errors due to time discretization and 
imperfect sampling could be diminished by the use of total sample sizes of 100,000 walkers and time steps of 0.001 
H-i. 



V. DISCUSSION AND CONCLUSIONS 



It should be mentioned that the relation between fermi statistics and the elimination of paths at permutation 
hyperplanes has been put forward by Korzeniowski et al. [p^ . This approach was subsequently criticized because 
it left open a few fundamental questions. Studying the underlying mathematical Ansatz has lead us to an extension 
of the recently reported many-body diffusion approach |^-^ in which it was the basic idea to rigorously separate the 
multi-dimensional free diffusion process of indistinguishable particles into one-dimensional free diffusion processes. 
The novelty of the outlined symmetry analysis for the simulation of indistinguishable particles is the construction of 
numerically implementable diffusion process with a switching of boundary conditions driven by a compound Markov 
chain . The continuous-time Markov chain, not considered before in that framework, has been set up to enable the 
realization of interdependent Brownian motions. The inclusion of transitions between a set of simultaneously evolving 
diffusion processes has been shown necessary to approach identical particles in potentials which are not invariant 
under permutations of Cartesian coordinates. To transfer the recently reported MBDA sampling technique [Q to those 
systems, a symmetry analysis of the potential is necessary. In applications to different model systems, the consequences 



14 



for the total diffusion process have been demonstrated. The numerical findings confirmed the feasibility of the utilized 
symmetry concept. The model of one particle exposed to a double-well or a Mexican-Hat potential transparently 
illustrates the essence of our approach. Artificially creating a sign problem by the decomposition of the propagator in 
two resp. four subpropagators, some of which are not positive definite on the configuration space, the efficiency of the 
algorithm has been found comparable to the corresponding boundary-free formulation. This application points out 
the generality of the symmetry analysis on which the algorithm is founded; its use is not restricted to the description 
of identical particles but may also be of relevance for the implementation of boundary conditions on distinguishable 
particle systems. In the MBDA, the boundaries guarantee the generation of symmetric and antisymmetric Brownian 
motions. This is fundamentally different from nodal approaches; obviously, the ground-state wave function of the 
particle does not vanish at the boundaries (cfr. figs. 2 and 4)! The flexibility of the many-body diffusion Ansatz 
has been indicated by treating two identical fermion oscillators in three dimensions with an algorithm based on the 
same principles but different subprocess-transition probabilities. Evidence has been given that the MBDA avoids the 
fermion sign problem in applications to the ortho-helium model. Necessary ingredients for this study have been the 
derivation of the appropriate reduced state space with the properly adapted stochastic processes. The comparison 
of the numerical ground-state estimates indicated that the MBDA reliably simulates the exchange contributions. Over 
the considered evolution intervals no significant decay in the signal-to-noise ratio has been observed. The reduced 
state space boundary 9I?2 does definitely not represent a nodal surface. Also in numerical practice, the characteristic 
fermion state is completely determined by the Euclidean-time Schrodinger equation supplied with the symmetry 
properties of the free-diffusion process and the potential under consideration. From that point of view, the MBDA 
represents a potential candidate for ab-initio fermion Monte Carlo simulations. In comparison with literature the 
MBDA — in its present form — is restricted by its versatility. The MBDA tackles fundamental methodological aspects 
to avoid the sign problem for fermions. To transparently demonstrate that the supply of nodal surfaces is unnecessary 
for the current approach, the present formulation dispenses with such analytical assumptions. We are well aware 
that importance sampling may drastically enhance computational efficiency, but the objective has been to keep the 
algorithm as transparent as possible. The nodal structure of the ortho-helium ground state, for instance, can be 
derived by symmetry arguments [T^ . Implementing an importance function with this nodal structure enables one to 
circumvent the sign problem. The MBDA algorithm, however, does not require to build in the nodal structure. We 
stress that the MBDA is fundamentally different from approaches based on interacting or correlated walker ensembles 
p8|-p2|. Dealing with many- fermion diffusion on the appropriate state space Z?, the ground-state wave function is 
simulated by a sum of antisymmetric elementary densities. This is in contrast to the common definition of positive 
and negative walker ensembles, the densities of which behave naturally symmetric. The MBDA involves only positive 
walkers. Antisymmetry is realized by establishing the correct free diffusion process without requiring projection by 
an antisymmetric trial function. The use of reflection and absorption restricts the evolution to the state space; this 
procedure proves useful as long as the supplementary effort for reflection and absorption is balanced by enhanced 
accuracy due to restriction to the state space. The present investigations mainly attempted to solve and clarify the 
fundamental difhculties related to the sign problem for fermions. The underlying symmetry concepts were analyzed 
and used to guarantee the stable simulation of fermion systems. Some transparent examples were chosen mainly for 
illustrative purpose and no inference has been made to provide ultimate accuracy. 
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Figure Captions 



Fig. 1 : The rigorous (right hand side) and numericaUy sampled (left hand side) Green function of two free identical 
fcrmions in one dimension restricted to the reduced state space xi > \x2 \ after an evolution of three atomic time units. 
The initial condition is denoted by the crosses, the state-space boundaries are indicated by the solid lines. 

Fig. 2: Plot of the symmetric and asymmetric double well potential Vs{x) and Va{x) ( p7|) with a = 2. Solid 
line: 6 = 0, dashed line: b = 0.25. The long-dashed resp. dot-dashed curve depicts the corresponding simulated 
ground-state wave functions ^s(a^) and 'I'a(a;). 

Fig. 3: Plot of the Mexican Hat potential ( p9| ) with a — 2 and b — (l.h.s.) and with a — 2 and b — 0.25 (r.h.s.). 

Fig. 4: Comparison of the two numerically predicted ground-state wavefunctions of the asymmetric Mexican-Hat 
model (29), with a — 2 and b — 0.25 sampled with the MBDA (l.h.s) and the boundary-free algorithm. 

Fig. 5: The ground-state energy estimates (left hand side) and the relative walker population of the subdensities 
(right hand side) as a function of Euclidean evolution time for the three-dimensional anisotropic system of two identical 
fermion oscillators. 

Fig. 6: Comparison of the numerical energy estimates (solid curve) for the lowest ortho-heliimi state with the 
numerical estimate -2.175229 H found in ||l^ (dashed curve). The dot-dashed lines denote a deviation of one percent 
from the dashed curve. 

Tab. 1: Estimated ground-state energies Eq with standard deviations a and relative subprocess population for 
the double- well potential potential (27). 

Tab. 2: Same as table 1 but for the Mexican-Hat potential (pSI). 
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